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An approximate analytical and non-linear solution of the Einstein field equations is derived for 
a system of multiple non-rotating black holes. The associated space-time has the same asymptotic 
structure as the Brill-Lindquist initial data solution for multiple black holes. The system admits 
an Arnowitt-Deser-Misner (ADM) Hamiltonian that can particularly evolve the Brill-Lindquist so- 
lution over finite time intervals. The gravitational field of this model may properly be referred 
to as a skeleton approximate solution of the Einstein field equations. The approximation is based 
on a conformally flat truncation, which excludes gravitational radiation, as well as a removal of 
some additional gravitational field energy. After these two simplifications, only source terms pro- 
portional to Dirac delta distributions remain in the constraint equations. The skeleton Hamiltonian 
is exact in the test-body limit, it leads to the Einsteinian dynamics up to the first post-Newtonian 
approximation, and in the time-symmetric limit it gives the energy of the Brill-Lindquist solution 
exactly. The skeleton model for binary systems may be regarded as a kind of analytical counterpart 
to the numerical treatment of orbiting Misner-Lindquist binary black holes proposed by Gourgoul- 
hon, Grandclement, and Bonazzola, even if they actually treat the corotating case. Along circular 
orbits, the two-black-hole skeleton solution is quasi-stationary and it fulfills the important property 
of equality of Komar and ADM masses. Explicit calculations for the determination of the last stable 
circular orbit of the binary system are performed up to the tenth post-Newtonian order within the 
skeleton model. 

PACS numbers: 04.25.-g, 04.20.Fy, 04.25.Nx 



I. INTRODUCTION 

The description of the motion of binary black holes 
within the Einsteinian theory of gravity is quite a chal- 
lenging problem and is likely to be finally solvable by 
numerical means only. However, it seems that the full 
success will occur in a far future. It has turned out in 
the past that analytical developments did always influ- 
ence numerical investigations, hence further progress on 
the analytical side is certainly desirable. In this paper we 
shall present an analytical solution of truncated Einstein 
equations for systems of non-spinning point-like objects, 
which shows many aspects of past numerical models. 

Most of the numerical computations of quasi- 
stationary initial data for black-hole binaries are based on 
the assumption that the space metric is conformally flat 
USQ. They often assume as well the existence of an 
approximate helical Killing vector which permits to elab- 
orate some manageable formulations (also see 4J). Re- 
cently, the time-symmetric initial data of Brill-Lindquist 
P| were generalized by adding a non-conformally flat 
contribution that incorporates pieces of information pro- 
vided by perturbative post-Newtonian calculations 0. 
The present paper does not follow this line, but rather 
keeps the spirit of former numerical simulations. We es- 
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sentially perform two simplifications: on the one hand, 
we adopt the often-used assumption of conformally flat 
space metric; on the other hand, we attribute the non- 
linear gravitational field energy to point-like sources, 
which are modeled as Dirac delta distributions. We solve 
the resulting truncated Einstein equations analytically 
and get what we call, in reference to the point-like char- 
acter of the source support, the skeleton solution. The 
calculations are carried out for Af bodies in a mathe- 
matically sound way by working from the beginning in 
(e?+l)-dimensional space-time and taking the limit d — > 3 
in the final expressions, so that the Dirac delta distribu- 
tions modeling the sources are consistently handled 0- 

Our space metric and extrinsic curvature may be 
regarded as some generalization of the multiple black 
hole Brill-Lindquist initial value solution p| to arbitrary 
black-hole momenta preserving its original asymptotic 
structure. Moreover, the energy function can be used as 
generator of the dynamical evolution of the system. We 
shall refer to it as the skeleton Hamiltonian associated to 
the full Arnowitt-Deser-Misner (ADM) Hamiltonian. It 
agrees with the exact Hamiltonian (i) in the test body 
limit, (ii) at the first post-Newtonian (1PN) approxima- 
tion, i.e., at the first order in powers of 1/c 2 (c denoting 
the speed of light) beyond the Newtonian dynamics, and 
(in) in the limit of special relativity in absence of gravity. 
The skeleton gravitational field is given explicitly, as well 
as the equations of motion. 

As a first application, we shall investigate the loca- 
tion of the last stable (circular) orbit (LSO) for binary 
systems. In order to reach a proper accuracy, the last- 
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stable-orbit parameters are worked out up to the tenth 
post-Newtonian (10PN) approximation within the skele- 
ton model. The comparison of second and third post- 
Newtonian results derived from the skeleton model with 
the corresponding quantities computed in the full Ein- 
stein theory delivers a quantitative measure for the miss- 
ing terms in the skeleton approach. 



II. DESCRIPTION OF THE BLACK-HOLE 
DYNAMICS IN ADM FORMALISM 

We are interested in a system of Af black holes inter- 
acting gravitationally. Its description will be achieved 
within the ADM Hamiltonian framework using a 1) 
splitting of space-time. The spatial indices are denoted 
by Latin letters and vary from 1 to d, whereas the space- 
time indices are Greek and vary from to d. In a given 
chart, (x = x l , t = x°/c), the motion of the Ath black 
hole is described by a curve x A {t), for A = 1, . . . , N. 
The conjugate momentum of the position x\ is called 
PAi- The constant mass parameter of the Ath black hole 
is denoted by ttla- The Hamiltonian variables associated 
to the gravitational field are taken to be (i) the space met- 
ric jij induced by the full space-time metric g M „ on the 
hypersurface t = const, and (ii) its conjugate momentum 
(c 3 /WirG)^ where G is the Newtonian gravitational 
constant. The other components of t/ M „ are parametrized 
by the lapse and shift functions: N = {goigoj'f^ — goo)^ 
and Ni = gm, with 7 y standing for the inverse of jij. It 
is also useful to introduce the shift vector N l = 7 y iVj. 
The spatial indices will be raised and lowered with the 
help of r ) lJ and 7« , respectively. 

We now assume that the "center-of-field" x^(t) of ob- 
ject A can be modeled in the fictitious space-time asso- 
ciated with the flat metric by means of the Dirac distri- 
bution 5a = 5(x — x^). This prescription is known to 
lead to the Brill-Lindquist Hamiltonian when the black 
hole velocities are set to zero, provided the mass param- 
eters tua, A = 1, . . .. Af, are identified with the Brill- 
Lindquist masses 0, [jj ■ It also leads to the correct 2PN 
dynamics |Tol ITU for compact binaries, which can itself 
be recovered by means of the extended-body approach 
[l2|. Moreover, the singularities generated by the dis- 
tributions 5a have shown to be curable by dimensional 
regularization at the 3PN order 0- Similarly, we shall 
see in the present work that, for an appropriate choice of 
the space dimension d, the <5A-type sources do not gener- 
ate any ambiguity, so that the 3-dimcnsional quantities 
can be obtained by dimensional regularization. Finally, 
this point-particle approach allows to bypass delicate in- 
vestigations about the behavior of the gravitational field 
near the objects. 

In an asymptotically flat space-time of dimension d+1, 
the Hamiltonian of Af point masses provided by the ADM 
formalism reads [lj, Ej, EH 



H=^[jd^{NH + N^) 

+ J d"x i) Mir,, n,- u ^. (2.1) 

where di = d/dx l . The super-Hamiltonian density Ti. 
appearing under the integral sign is defined by 

j\f i 

16ttG ^ ( 2 ^PAiPAj \ 2 , ,„ v 

+ 2^[ m A + ~2 ) S A, (2.2a) 

while the supermomentum density Ji is given by 

16irG ^ 

Ji = -20^\ + 7T kl d ilkl — PmSa- (2.2b) 

c A=i 

In Eq. H2.2afl R denotes the space curvature of the hyper- 
surface t = const and 7 is the determinant of the matrix 
7ij . As the lapse and shift functions are Lagrange multi- 
pliers, the variation of the Hamiltonian (|2.1|) with respect 
to them leads to the constraint equations 

H = 0, Ji = 0. (2.3) 

Before solving Eqs. I|2.3I) . we must first fix the coor- 
dinate system. It is convenient to work in the so-called 
ADM transverse trace-free (ADMTT) gauge [H H3 de- 
fined by the two conditions 

dj (la - ~<%7fcfe) = 0, n u = 0. (2.4) 

In this gauge, the metric decomposes into its trace 7^, 
which may be parametrized as d\E» 4 /( d_2 ) (in order to 
get a simple expression for the curvature density), and a 
transverse trace- free part hj? . We thus have 

Jij • 9^=0, ^ T = 0, (2.5a) 

tt" = 0. (2.5b) 

The momentum n lJ needs itself to be split into a trans- 
verse trace-free contribution ir^ T , and a rest TT tJ , 

n ij =ff ii +7r « r _ (2.6) 

These two terms are uniquely defined assuming they 
decay as 1 / (x l x l )~%~ at spatial infinity and are given 
as a combination of derivatives of the Poisson integral 
A~ 1 djir l: > , and the Poisson integral of the latter quan- 
tity, namely A~ 2 djir l: > . 

We can then try to compute and ?r y from the two 
constraint equations (|2.3|l supplemented by the boundary 
conditions at spatial infinity resulting from the asymp- 
totic flatness of space-time; notably \& — l^(x k x k )^~ , 
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K l i ^l/(x k x k )^~ . It is actually more convenient to 
solve the Hamiltonian constraint for the quantity <f> = 
^tEz — -0 rather than $ to get rid off the spurious 
constant at spatial infinity. We thus have 



* = 1 



d- 2 
4(d = l)' 



After solving the constraint equations (|2.3|) . we insert 
the ensuing expressions for (f> and into the right-hand 
side of Eq. (|2.1I) . We get a reduced Hamiltonian depend- 
ing on x\, pAi 



hj? , and 7r^ T only: 



H lc d — i? rc d[s/i, PAi j hij^i ^tt] 



16nG 



(2.7) 



This Hamiltonian contains the full information about the 
evolution of the matter (x l A ,pAi) and field (/i^ T ,7r^ T ) 
variables 0, . As the change of variables (7^ , ) — > 
(hj? ,tttt) is n °t canonical, the Poisson brackets of /i^ T 
and 7r^ T take a form involving a transverse trace-free 
projection operator 0, . 



III. SKELETON HAMILTONIAN 



hJT to be zero on £0 is always possible, though it is not 
physically realistic in general for isolated self-gravitating 
systems. Our attitude in this instance will be to accept 
the resulting error. In return, we shall be able to deter- 
mine the field 7y without assuming zero velocities of the 
bodies at initial time. 

We thus demand here that 



hj? = 0, 



(3.1) 



and the problem reduces to solve the constraints 7i = 0, 
Ji = under this condition. The space metric 7^ is then 
proportional to &y, which implies that in the supermo- 
mentum J7i the term w^d^ki vanishes [cf. Eq. Ij2.2b|) ]. 
The field momentum 7r y relates to tt*. through the iden- 
tities 



(3.2) 



In particular, 7r ! ? - is symmetric and trace-free (STF) as is 



The determinant 7 appearing in 7i is merely 

'y = \]/ d-2 



(3.3) 



Finally, the 3-curvature density R^/7 takes a remarkably 
simple form |22^ : 



= -$A4>. 



(3.4) 



Collecting these relations together with the constraint 
equations (|2.3|) . we find: 



The reduced Hamiltonian H rc a provides the equations 
of motion of the Af bodies together with the evolution 
equations of the gravitational field. As they are strongly 
coupled it is hopeless to solve them analytically, but they 
can be solved perturbatively us ing , notably, the post- 
Newtonian iterative scheme ^3 0, ^| . The resulting 
dynamics is highly reliable as long as the typical speeds 
remain smaller than, say, about 0.3c. When binary black 
holes reach their last stable orbits, the reliability of the 
approximation breaks down since the system enters a 
strong field regime. The only way to understand precisely 
this crucial stage is to perform numerical simulations in 
full general relativity. Our purpose here is to search for 
a simple non-perturbative model applicable to arbitrary 
strong fields, eventually leaving the Einsteinian theory 
but always staying close enough to it, in order to give ^ (part) 
an element of comparison with numerical investigations. 
The space metric obtained in our approach should also 
be usable as initial condition. 

We shall adopt here a Wilson-Mathews-type prescrip- 
tion [2(], , setting the non-conformally flat part of 7^ , 
namely hj? , equal to zero. As hj^ is of order 1/c 4 , the 
difference between our 7^' and the Einsteinian space met- 
ric shows up at the 2PN level. For spherically symmetric 
space-times the difference is zero. In the case where the 
present model is used to fix initial Cauchy data, the con- AVi = 
dition of a conformally flat space metric needs to hold 
on only one hypersurface, say Eq : t = to- Choosing 



6 A=l 



ttiaSa^' 



t-~£- Pa 

V]> d-2 A 



d,K = -—^-2^PAi5 A , 

C A=l 



(3.5a) 
(3.5b) 



With p\ = PAiPAi- 

In the system lj3.5[l the second equation decouples from 
the first one and is solved by searching for a particular 
solution proportional to the symmetric trace-free part of 
the derivative of a certain vector potential Vi] namely, 



STF(2d i V j ) = d i V j +d j V i -~6 ij d k V k . (3.6) 

The equation for the divergence dkVt is obtained by ap- 
plying the operator dij = didj simultaneously on the 
two sides of (|3.6|1 . hence 2AdkVk = d/(d — l)9y7r^. 

This yields immediately djir{ = AVi + (f - l)/(d - 
l)A~ 1 dijkTr 3 k , and from the constraint (|3.5b(l it follows 
that: 



8ttG / d — 2 1 

■—^[pmSa-^j^paAA ^ 



(3.7) 
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We now notice that the <i-dimensional Laplacian of the 
mth power of = [{x l — x l A ){x l — x^)] 2 satisfies 



Ar r ; 



m(m — 2 + d)r T j 



(3.8) 



in the sense of functions, and that, in the sense of distri- 
butions, 



2-d 



Ar 



As a result we may choose: 



47T2 



-6 A - 



(3.9) 



Vi 



2tt- 



4gAj 
„<2-2 



(d - 2)pAj 
(d-l)(4-d) 



»4-d 



l J A 



(3.10) 

which is nothing but the 1PN vector potential V? 3 n of Ref. 
verifying n ij = STF(2diV ( L) + C(l/c 5 ) 



Our second simplification consists in stating that the 

m J- ( h ° m ) — 



homogeneous solution of Eq. 



(part) 



is precisely zero, which implies 7r 



i TT 

3 



0. 



When this prescription is imposed on one unique hyper- 
surface, it merely represents a particular choice of initial 
data. The resulting field momentum = STF(29jV}) is 
identical to the rescaled extrinsic curvature proposed by 
Bowen and York in their time-asymmetric formulation of 
the Cauchy problem |2^| and adopted e.g. in the punc- 
ture method l24| or various calculations of the last stable 
circular orbit [2J, |2o> |2]| . Although it is a mathemati- 
cal solution of the constraint equations for a conformally 
flat metric, it is not realistic in the sense that it does not 
result from an evolution process in absence of incoming 
radiation. It does not coincide either with the Wilson- 
Mathews momentum (7t 1 )wm derived from the condition 



L TT 



= 0: 



(ttI-)wm 



(1- 



i (part) 



, (3.11) 



where the quantity defined as the 1/c 5 part of 

the solution [matching the form (3.6)] of the equation 
(d — l)djTt%s = — dj(</)7r J i < ' part ^), is explicitly given in pa- 
per 0| for d = 3. The difference (t^wm — 7Tj ( part ) i s 
only of order 1/c 5 , which corresponds to the 2PN approx- 
imation. On the Hamiltonian level, it is of 3PN order 
only. One might then argue that the Wilson-Mathews 
condition hjj 1 — 0, valid for all times, i.e. particularly 

hjj 1 = 0, is in contradiction to the condition 7r^ TT = 0. 
However, it is not the case because the momentum con- 
straint equation is linear in Trt and decouples from any 
other field equation when the space metric is conformally 
flat, so that the quantity obtained by adding an arbitrary 
homogeneous solution to tt 1 - still satisfies the constraint. 
As our it 1 - differs from (tt*-)wm by such a homogeneous 
solution (i.e. a transverse-trace free object), there is no 
contradiction between the two approaches. 



The integration of the Hamiltonian constraint (|3.5a() 
is much more involved, principally because of the ap- 
pearance of <p in the source and the presence of a non- 
compact support term. On the contrary, the compact 
support term is harmless. It is made of a sum of Dirac 
distributions 5a that formally determine the boundary 
conditions near the black holes. By subtracting from the 
conformal factor the 1/r^ -2 terms entering (/>, we get an 
object that satisfies an equation whose source is entirely 
free of Sa's (this is the essence of paper j^). By con- 
trast, we shall keep here the constraint equation (|3.5a(l 
under its original form in order to determine the depen- 
dence of the coefficient of l/r^~ 2 on the positions x\ and 
momenta pAi- 

To solve the Hamiltonian constraint 1(3. 5aj) itera- 
tively, we have to face the delicate computation of 

A I W d ~ 2 7r*-7Tj J . The expansion of the prefactor 



2-3d 

qj— = 1 1 



2 

d-2 



in powers of <fr is known to 

give rise to poles when d — > 3 at the level correspond- 
ing to the 3PN approximation 0. In full general rela- 
tivity, these poles cancel in the Hamiltonian with simi- 
lar quantities coming from hj^ , but such a cancellation 
cannot occur in our case. In particular, the 3PN ap- 
proximation of the Wilson-Mathews model does not ex- 
ist for binary black-holes. We must therefore treat the 
7r 2 term differently. Let us first note that this term 
enters the space metric at the 3PN order only, since 
7^-7^ = C(l/c 6 ), but contributes to the 1PN Hamil- 
tonian [because of the global factor c 4 /(167rG) in Eq. 
(12.7(1 ], To be consistent with the accuracy of the for- 
mer assumption hj^ — which amounts to neglect some 
2PN corrections in the Hamiltonian, we have to keep 
at least the leading term in the post-Newtonian expan- 

2-3d 

sion of ^1 + 4 ( d rf l 2 i) <frj d 2 7r j" 7r i' modulo a possible total 

space derivative. Now, the tensor density 71"^ is itself pro- 
portional to the symmetric trace-free part of the space 
derivative diVj, so that 



2STF(7r l ,)STF(d 4 V,-; 



2n)diVj . 



As a consequence, the latter expression can be split into a 
"skeleton" term — 2VjdiTT l j with compact support, and a 
"flesh" term 9j(2Vj7Tj) which takes into account the field 
between the point particles. The flesh term contributes 
to the Hamiltonian through the quantity: 



16ttG 



d d x* — d l (2V j n)) 



3d -2 



32nG d - 1 



4(l-d) - ' 

*~* = *~Vjir)di<l> 



We are allowed to regard it as negligible in our scheme 
and keep the skeleton term alone in the decomposition of 
7r*-7r^. Our third simplification consits thus in performing 
the substitution 
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T i=M i i 
■qr^=TTr.Tr- 



-2*" 



= 1 -^^J2 PAjVj 8 A . (3.12) 

° A=\ 

After the replacement Ij3.12|) is made, the right-hand 
side of Eq. (|3.5a|l takes the form J2 A f A (x)5 A , where the 
f A s are some unknown functions. For d belonging to 
an appropriate range of values, f A is regular at x = xa, 
and A(j) = ^2 A f A (x. A )5 A . We see that in our model the 
source of the Hamiltonian constraint reduces to some 
"skeleton" made of a linear combination of Dirac dis- 
tributions. For this reason, all quantities computed in 
the present approximation will be referred to as skeleton 
quantities. Defining the Af quantities a A by the relation 



16ttG 



oa = Ia^a), A = l, 



the function <b can be written as 



4GT(^)^ aA 
2^ r d-: 

A=l A 



c Tt^r- 



(3.13) 



and the conformal factor reads 



OA 



A=l r A 



(3.14) 



where we have posed, in order to lighten the notation, 
Gd-2T(^) 



A,, 



c 2 d-l n ^ 



It is also useful to introduce the (regularized) value of W 
at the particle position x^: 



* A = *(x = x A ) = 1 + A d — 



a B 



B^A 



' AB 



where r AB = [(x A - x B )(x\ - x l B )\ 5 . 

After making the replacement 13.12fl on the right-hand 
side of Eq. (|3.5a|) , inserting there expressions (|3.13l) and 
(|3.14l) . we arrive at 



16ttG 



A=l 



Acj)=J2 vaSa = { m A ( 1 + A d -IT 



B^A T AB 



i+(i + A d y, 



OLC 



C^A ' AC 



a B \ d ' 2 p Ai V Ai 



.d-2 



>5 A , 



B^A ' AB 

with Vm = Vi(x = x A ). Making use of Eq. I|3.10[) . we see that 



PAiVAi 



A 



/3d -2 



ByiA V 



(PAiPBi) + (d-2){n AB pAi)(n AB pBj) )r 



„2-d 
AB ■ 



' ft 

9 9 

m A c z 



(3.15) 



(3.16) 



Because of the linear independence of the distributions 5a, the Hamiltonian constraint l|3.15(l is finally equivalent to 
a system of Af algebraic equations which reads 



ola = 



m A 



1 + Ad J2b^a r d-a 



P A /( m A c ) 



1 + A d J2c^A r "-2 



PAiVAi I C 



1 + A d ^2 B _ tA ^ 



3d-2 5 
d-2 



A = l, 



A/". (3.17) 



Solving the above system provides in principle the values 
of all the a A s. These quantities represent some effective 
masses, as suggested by the approximate equality a A = 
m A + 0(l/c 2 ). They are functions of x\ — x\ (for A ^ 1 ) 
and p A i [defined implicitly by Eqs. (|3.17|) ] . 
The skeleton Hamiltonian reads 

r 4 r M 

H = -— — / d rf x A0-c 2 Va A , (3.18) 



which shows that the ADM mass Madm = H/c 2 is just 
the sum of the Af effective masses a A . The Hamilto- 
nian l|3.18(l depends only on the matter variables: H = 
H(x l A — x\,p A i), so that the system is automatically con- 
servative. The space metric is deduced from Eqs. (|2.5all . 
(I3~T1) . and J3~T4l : 

\ A =i r A ) 
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By construction, H and 7^- coincide with their coun- 
terparts in the Einsteinian theory up to the 1PN order. 
The test-body limit is achieved by taking the mass of one 
black hole, say of label 1, much larger than the masses of 
its companion: mi 3> tua, where A = 2, . . . , TV. To 
get the test-body Hamiltonian exactly, one has to let 
formally mi go to infinity (mi — ■* +00) in the Hamil- 
tonian (|3.18() while keeping mi/r", A — 2,...,Af, 
constant as well as all linear momenta pAi- It follows 
that particularly the flesh term 9j(2Vj7T*-) does not con- 
tribute to this approximation. In absence of gravity 
(G — > 0) our Hamitlonian has the form of the exact 
special-relativistic Hamitlonian for the system of N free 
point particles. When pAi = for all A, we recover the 
Brill-Lindquist solution, for which Eqs. (|3.17ll reduce to 

aA(l + Ad ^2b^a a B/ rd AB) — m A- As a consequence our 
model also generalizes the Brill-Lindquist data. 

Let us observe that the particle masses tua cannot be 
equated to the ADM masses computed at the spatial in- 
finities of the AT black hole sheets when one or several 
momenta are non-zero. When such an identification is 
performed, it leads automatically to the same relation be- 
tween a a and tua as in the zero momentum case, whereas 
the correct dynamics up to the 1PN order is known to 
be the one that derives from the Hamiltonian 13.1811 de- 
termined by the set of Eqs. (|3.17l) with pAi 7^ 0. The 
actual dependence of the effective masses on the to^'s 
must take the kinetic energy of the system into account. 
This indicates a relative boost of the black hole sheets 
with respect to each others. The issue is complicated by 
the fact that the Bowen-York data for one single black 
hole do not result from a boosted Schwarzschild field [2^] . 
This point, being still rather unclear, needs to be inves- 
tigated in detail. We do not intend to discuss it here. 

1 We use here and below the usual notation for symmetrization: 
2D {i N j} = /.». .V, • D,S;. 



IV. SKELETON GRAVITATIONAL FIELD 

In the ADM formalism, the evolution of the space met- 
ric in (d+ 1) dimensions is given by [cf. Eq. (4.70) in Ref. 

where the operator denotes the space covariant deriva- 
tive. Its explicit action on the shift function after 
symmetrization is 

D(iNj) = lHi d j} N k + \N k d kll3 . (4.2) 

4 

For the conformally flat 3-metric 7y = ^> d - 2 5ij, we 
have 7y - \~ikk&ij = 0, hence 

d t (lij ~ -TkkSijj = 0. 

After replacing in the above equation dtjij and dtjkk by 
the expressions following from Eq. i|4.1[l and dropping all 
contributions proportional to ir k k — 0, we arrive at: 

STF(diN j ) = -V-^Ntt) . (4.3) 

The evolution equation for the field momenta ir lJ is 
known in d dimensions as well [cf. Eq. (4.71) in Ref. |l5|]. 
It reads 



-d t ir 13 = - V7 



N ( R y - - 7 y R ) - D l D J N + ^D m D m N 



N 



\^D k N* + n ki D k W ~ D k (n^N k )] + ^-N Y £4*Pd£ 7 <*y« f l + 



c — mA 



PAmPAn 

? 9 ' 
m A C Z 



6a- 



(4.4) 



In the ADM gauge, the trace (with respect to flat metric Sy) of 14.411 reduces to 



d- 2 



i-d 



NR — (d— 1)$-*^ DiDiN ) - ^^y^^ Ntt)-k{ - 2n kl D k N l 

87rG , T T _ s p 2 A 1 , 
-I — Af 3=2 \ ZA. 1 

c 4 f-^ m A » 

A=l 



2 _2 



m A c 



6 A = 0. (4.5) 



Making now use of Eqs. (|3.4|l . (|3.5a|) . (|4.3H . and of the two useful relations D^D^N = j k ^dj)(D k N) + ^d k j lj D k N 
and DiDiN = 2d i Nd l In * + AN, Eq. JOI can be rewritten as 
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4 c \ m A c ) 

+ ^-n N jr A( 1+ J&*-AY i s A = 0. (4.6) 



A=\ 



The form of the latter equation suggests to solve it for the auxiliary function \ = rather than for the lapse N |(| . 
At this stage, we apply to the right-hand side of Eq. (|4.6|) the replacement <|3.12|) which amounts to deleting a flesh 

3d-2 . ■ 

term in the product \& d ~ 2 ir^n -\. Since the Newtonian potential comes from the 1/c part of the lapse, the term 

^ -1 $ _ "3=3 _ 7rV7rj = C(l/c 6 ) is indeed a 2PN quantity and may be neglected consistently with our approximation. 
Finally, we obtain a Poisson equation for \ with the source term being some linear combination of the Dirac deltas: 



A=l \ 



d-2 

A 



Following the same argument as in the derivation of the 
conformal factor in Sec. IIIII we may write 



(4.8) 



so that Ax = ^gr- Ea=i We next replace the 

left-hand side of Eq. (|4.7() by virtue of the latter relation, 
and equate the coefficients of 8 a- We are led to 



Pa = 



B^A ' AB 



m A *~ A 2 [ 1 



m A c 2 A 



d-2 c 



( 1+ /„ P 2 \ ^ dl2 
d-2 m A c 2 A 



A=l,...,Af. 

(4.9) 



This system of J\f equations for J\f unknowns provides 
in principle the values of the monopoles Pa entering the 
expression of x, Eq. 1)4. 8|) . The lapse function itself is 
then given by 



For M = 2 and d — 3, we recover the result of Ref. [3 in 
the time-symmetric limit pAi = (with A = 1,2). 

In order to solve Eq. (|4.3|l for the shift function N % , we 
take the first and second order space derivatives of both 
sides. We obtain: 



„ / 2d 

-2dj — 



~—N-k'\ 



(4.10a) 



-A(a,iV J ) = -S!y (*-^JVtt*-) . (4.10b) 



The source term of the second equation is the 
sum of a rational function of ta , A — 1 , . . . , J\f, 
namely — di(dj( i S~' 3=5 N)ir l j) , plus a distribution, 

—di^^^^NdjH 1 ^. The Poisson integral of the term 
with non-compact support is not explicitly known, but 
that of the distributional term is easy to compute. Now, 
the new "flesh" term does not contribute at the Newto- 
nian level. It is therefore at least post-Newtonian and will 
eventually play a role in the dynamics only at the 2PN or- 
der; it will be consistently neglected in the present model. 
Thus, the substitution l|3.12jl will be supplemented by the 
replacement: 



dj Nitty -> V-^NdjTc'i 

8nG 



AT 



N^ PAi 6A. (4.11) 



A=\ 

Taking Eq. I|4.11|) into account, Eq. (|4.10b|) becomes 

Q_/^ J / N 3d~2 

A(d^) = ^ — dj £ XA* A d " PmSa 

C ^ A=l 

(4.12) 

The solution of 14.12J1 can be obtained with the help of 
formula (13.911. It reads 



2G d r(^) 
c 3 d - 1 ^ 

" _Mz? / 1 

£ XA* A " 2 PAidi — 
A=l VA 



(4.13) 
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The Poisson equation (|4.10a|l after inserting the above ex- 
pression for djN 3 and performing the substitution (|4.11|) 
takes the form 

i _2G d-2T{^) 

C 3 d - 1 



x 



XA^A^PAjd, 



A=l 



16ttG 



AT 



^2 XA* 



3d -2 
d-2 



PAi$A ■ 



A=l 



(4.14) 



By virtue of the formulas (|3.8|l and (|3.9() , the solution of 
Eq. 14.14JI is simply 



where the skeleton Hamiltonian H expressed in terms of 
the quantities ^a equals 



H = 



2c 4 r 12 
G 



(*i + *a-2). 



(5.4) 



Using the theorem of implicit differentiation, we put 
Eqs. (|5.3|) in a form apparently more convenient for nu- 
merical integration. We start by differentiating Eq. 1)5. 2|) 
with respect to the positions. We find (here B — 1,2 and 
A + B): 



d 2l ^ B = 



dx\ 



Vm ~ CaOu^a, (5.5a) 
VA2 ~ (Ad 2 i*A, (5.5b) 



r 3 §=1 



A=l 



d-2 



4-d 



(4-d)(d-l) PAj ° ij1 " A " " ^ 



(4.15) 



V. BLACK-HOLE BINARY EQUATIONS OF 
MOTION 



For a black-hole binary, the number of coupled equa- 
tions composing the system l|3.17[l that gives the effective 
masses reduces to two, and the sum J2c^A a c l r AC con " 
tains only one term. It is now convenient to regard the 
two quantities 'J ,4 as the unknowns: 



OL B 



C^A r AC 



A =1,2, B^A. (5.1) 



Let us restrict ourselves to the space dimension d = 3. 
The constant Ad reads then A 3 = G/(2c 2 ), and the sys- 
tem <|3.17ll consists of two equations for the two unknowns 



where 



Vai 



_ GmAn\ 2 



VA2 

(a 



Pa 



2c 2 r\^ A 
Gdii{pAjV A] /ri 2 ) 



m A c 2 ^ A 



2c^\ 



Gm A l + 3 P 2 A /(m A c 2 ^ 4 A ) 
'2c 2 r 12 * 2 A [l+p A /(m A <?* A p 

7G pAj Vaj 
2c 3 r 12 ^ A 

Xb - 1 
XA 



(5.6a) 
(5.6b) 



(5.6c) 



Equality (|5.6bf) follows from the fact that the right-hand 
side of relation (|5.2|l depends on x\ and x\ only through 
the black- hole separation r 12 = [(x\ — x\){x\ — x 2 ))^ , and 
that we have dur\ 2 = —d 2 ir\ 2 . Solving Eqs. (|5.5a|) with 
respect to du^B and Eqs. (|5.5b|l with respect to Bi 
we find 



*r = 1 



GniA 
2r 12 c 2 ^ A 



PA 



A^A 



GpAiVAi 

2r 12 c^\ ' 



5 = 1,2, A^B. (5.2) 

It is clearly possible to decouple the two equations 1)5. 2J) . 
In Appendix ^] we show that each $b is a root of a 
polynomial of order 200. 

The equations of motion of a black-hole binary read 



X A 



dH 

dpAi 
dH 



(5.3a) 



1-CiCa 



(5.7a) 
(5.7b) 



It remains to evaluate the momentum derivatives 
d pu V B = d$ B /dpu and d P2 ^ B = d^ B /dp 2i . The sys- 
tem of equations for d Pli ^B is given by 

d Pu ^B = A i-Ud Pu ^A, B = l,2, A^B, (5.8) 
where we have posed: 

„i G PuSai 



Ml 



2m A c 4 r 12 ^ 



r,2 
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G d pii (p A jVAj) 
2c 3 r 12 ^ 7 A 



(5.9) 



The coefficients (a are identical to those appearing in 
dii^B because they represent the result of the differenti- 
ation of the right-hand side of relation (|5.2|1 with respect 
to ^>a at x\ and pAi constant. The solution of Eqs. (|5.8|) 
with respect to d Pli ^ b reads 



'At 



CaQ 1 bi 



1 - C1C2 



(5.10) 



from which we deduce d P2i ^ b by exchanging the particle 
labels. 

It is convenient to introduce the rescaled variables 



Va 



2c 2 ri 2 



G 



' Vaai 



>ab 



2c 4 n 2 ffi 
G 



' AB^ 



A = 1,2. 



(5.11) 



Making use of Eqs. (|5~3|) . (|5~7f) . lf5~TU|) . and (|3~TT|) . the 

equations of motion (|5.3|) can be written in the form 

i\ = T3^(( 1 " Ca)^! + (1 ~ Ci)^i) , (5.12a) 

pii = - r^-^ 1 - M - (1 - Ci)r?i) 



2c^ 
G 



+ * 2 - 2)ni 



2 

= r^^ 1 - - ( x - 

-c(^7 + ^)^«(py^i). ( 5 - 12b ) 

where the last equality follows from relations (|5.2|l and 
definitions (|5.6a|) - (|5.6b|) . Similar equations for x\ andp2i 
hold with the role of labels 1 and 2 exchanged. Knowing 
the positions and momenta at a given time t, we can 
solve the algebraic equations lj5.2|l numerically, calculate 
rf A , 9\ B and (a by means of Eqs. I|5.6[) . (|5.9|) . I|5.11[) . 
and determine x\ (t + dt) as well as pAi {t + dt) from the 
evolution equations (|5.12|l . 



adiabatically in time while the frequency increases ac- 
cording to the relativistic version of the Kepler law. The 
process ends as soon as the first instable orbit is reached 
before the coalescence phase. The point where this oc- 
curs changes according to the model (see Refs. [2^, |3(j 
and references therein for some comparisons). 

For binary black holes the skeleton Hamiltonian reads 
H = (ai + 0.2) (? 1 where ot\ and a 2 are solutions of the 
coupled system of equations l|3.17|l . It is impossible to 
solve this system analytically, but it is not difficult to 
solve it perturbatively within the post-Newtonian set- 
ting, treating e = 1/c 2 as a small parameter. At the nPN 
approximation, we search the post-Newtonian expansion 
of the coefficients a a in the form 



a a = m A 



n+l 



0(e 



n+2\ 



A=l,2. (6.1) 



Note that we need to carry the series up to the order 
7; + 1 in e to get an nPN-accurate dynamics with respect 
to the exact skeleton evolution (which first differs from 
the exact Einstein's evolution at the 2PN order). We 
now substitute Ijtj.lfl into Eqs. I|3.17fl and expand their 
right-hand sides up to the order n + l. By identifying 
the coefficients of different powers e l , i = 1, . .. , n + 1, 
we obtain a system of algebraic equations for aj^ . We 
display here the first three equations: 



(0) 
a i 

mi 



J 1 ) 



mi 



PuPii 
2m\ 



Gm 2 
2n? 



(PuPu) 
8mf 



Gm% 
2ri7 



OPliPli 

2ml 



2m\m 2 



n l 2 n \ 2 PuP2j a 2 



(o)' 



2m\m 2 



m 2 



Ar 2 



(2) 

mi 



(giiPij) 
16mf 



Gm 2 1 9(piiPu) 2 
2r 12 { 8mj 



r (0) 

2mlm 2 



m 2 



G 2 m\ ( Vop\iP\i 49p u p 2i 



z ' 12 



Ami 



Am\m 2 



VI. LAST STABLE CIRCULAR ORBIT 

As an example of application of our skeleton Hamil- 
tonian, we shall discuss in this section the parameters 
of the last stable (circular) orbit (LSO) in the relative 
motion of two inspiraling compact objects. This class of 
systems is of primary importance, being one of the most 
promising sources of detectable gravitational waves for 
ground-interferometry experiments as LIGO (Laser In- 
terferometer Gravitational Wave Observatory), VIRGO, 
GEO600, and TAMA300. Because of gravitational-wave 
emission, the separation between the bodies decreases 



ln\ 2 n\ 2 p u p 2j 
Am\m 2 



„(°) ' 



m 2 



G 3 m 



°'12 



similar equalities hold with labels 1 and 2 exchanged. It 
is straightforward to solve up this system to the required 
order. 

For convenience, we shall consider instead of H the 
dimensionless Hamiltonian (after ruling out the constant 
rest-mass contribution) 



H 



H — (mi + m 2 )c 2 
ac 2 
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fi — m,\vnil(m\ + rri2) denoting the reduced mass of the 
system. At the nPN level, we have 

2 n+1 

H = + O (e n+2 ) • (6.2) 



A=l i=l 



We restrict our study to circular orbits in the center- 
of-mass reference frame. Within the ADM gauge, the 
frame shift is achieved by imposing the condition 



Pli + Vii = , 
whereas circularity requires that 

n\ 2 p u = n\ 2 p 2i = . 



(6.3) 
(6.4) 



The form of H becomes particularly simple when use is 
made of the reduced relative position and momentum: 



G(mi + ni2) 



P, 



Pu 



P2i 



When both relations (|6.3() and 16.4fl hold, the Hamil- 
tonian (|(j.2l) can be expressed with the help of r l and 
Pi only. Notice also that it depends on the masses m\ 
and rri2 exclusively through the symmetric mass ratio 
v = mim 2 /(mi + m 2 ) 2 . 

In the next step, we introduce the conserved (reduced) 
angular momentum j along the cirular orbit of radius r: 



P 2 = — 



(6.5) 



After eliminating the linear momentum pi by means of 
Eq. (|6.5|l . the dimensionless energy E = H(r,j) of the 
binary system becomes a function of the variables r and 
j: 

E j2 1 j. f 1 u ji ( 1 .x. ^ V 

7-^--+U2-^( 3 + -) + i^(- 1 + 3 -))2 



2^3 (! + ") + ^(5 + + -4". 



(6.6) 



The coordinate radius r of the circular orbit and the 
angular momentum j are related through one of the 
canonical equations of motion, namely 



dH(r,j) 
dr 



0. 



(6.7) 



a relation that may be regarded as expressing the quasi- 
stationarity of the configuration. Remarkably, it can be 
shown in this case, i.e., for circular motion, that the 
ADM mass of the system equates the so-called Komar 
mass Mk, i.e., the mass- like quantity proportional to the 
monopole part of goo at spatial infinity satisfying the re- 
lation 5oo = -1 + 4A 3 M K /r + 0(1/ r 2 ) H3- The proof 



of the equality Madm = Mk is given in Appendix for 
an arbitrary number of space dimensions. 

Equation (|(j.7|) is solved perturbatively for r = r(j). 
By plugging the result into Eq. I|6.6|l . we find the gauge 
invariant link between the energy E and the angular mo- 
mentum j in the case of circular motion: 

E _ H{r(j),j) 

e e 



1 I 

1 



2 fi , 4 ( 9 + -) 



16/ 



(81 + 41i/-5i/ 2 )e 2 + C(e 3 ). (6. 



It is convenient to rewrite formula l|6.8|) in terms of the 
dimensionless variable 



dE 
If 



(6.9) 



The parameter x is now treated as a small quantity of the 
order e [cf. Eq. !|(j.8[) ]. The definition (|6.9I) establishes 
the link between x and j. We invert it order by order 
to obtain j as a function of x. Inserting j = j(x) into 
Eq. (|6.8|l leads to a gauge invariant expression of the 
energy E depending on the parameter x. The calculation 
is performed at the 10PN approximation in order to come 
reasonably close to the exact result. We display here the 
explicit formula with the 3PN-accuracy. It reads 



E 



x /3 v 
~ 2 + \8 + 24 

/675 8585 
+ \128 + "38? 

n 

i=S 



27 
16 



29 
— i 
16 



17 

48" 



v 2 )x 3 



7985 
192 1 



1115 

10368" 



^U 4 



(6.10) 



The higher order post-Newtonian coefficients e, (for i — 
5, . . . , 11) can be found in Appendix IO 

We want to compare the result of Eq. I|6.10|l with the 
relation Ewm = %m (x) of the Wilson-Mathews model 
(remember that in this model the post-Newtonian expan- 
sion terminates at the 2PN order) 2 : 



X 

2 



x 



27 
16 



1 v 
B + 24 
39 17 



16 



-v 



x 3 + 0(x 4 



(6.11) 



The 2PN coefficient proportional to v 2 differs from that 
of the post-Newtonian computations carried out in gen- 
eral relativity. Indeed, the energy for circular orbits in 



2 Thi s relation is obtained by taking n = j in Eq. (C4) of Ref. 
|32|l and replacing there the angular momentum j by jwM(%) 
according to Eq. (Bl) of Ref. |29| . 



11 



Einsteinian theory, up to the 3PN order, is [2 
'3 




27 19 1 

v H v" ]x° 

16 16 48 ' 



0{x% 



(6.12) 



where, according to Ref. we set the parameter w sta tic 
(defined in paper [3j|) equal to zero (w sta tic = 0), getting 
the solution corresponding to the Brill-Lindquist initial 
data. 3 By comparing Eqs. (|6.1U|) and l|6.12|l it is obvi- 
ous that our energy E starts to differ from the general- 
relativistic one, Eqr, at the 2PN order. 

The position xlso of the last stable orbit is deter- 
mined from the requirement that (dE / dx) I=ILSO = 0. 4 
The corresponding energy Slso [computed by means 
of Eq. H6.10fl ] is depicted in Fig. ^ f° r equal-mass bi- 
naries (i.e., for v = 1/4) as a function of the dimen- 

sionless orbital frequency Xlso- Numerical values are 
given in Table [fl The case of a test particle moving in 
a Schwarzschild background is also plotted for compar- 
ison. Similarly to what happens in Einsteinian theory 
[30|. the instabilities arise "earlier" (i.e., at lower fre- 
quency) for higher post-Newtonian orders in the skeleton 
model. The 2PN last stable orbit of our skeleton dynam- 
ics is quite far from that of general relativity in the sense 
that the absolute difference <5(^lso) 2PN W-^lso^pn) 
between the skeleton frequency (a^go^PN = 0.1077 
(the skeleton energy (-Elso)2PN = —0.0689) and the 



Einsteinian one (x%^ )2pn 

LSO 

—0.0795) is not negligible with respect to the general- 
relativistic value: Six^so) 2 ™/^^ )2PN = 21% 

LSO 

(6(E lso )2Pn/(E gr ) 2 pn = 13%). The location of the 
Wilson-Mathews data suggests that modifications of the 
extrinsic curvature affect the parameters of the last sta- 
ble orbit more than the choice of a conformally flat space 
metric ^3j. The 3PN skeleton values are even further 
away from the Einsteinian 3PN ones, but they are closer 
to those of the effective one-body approach improved by 
Pade approximants than the straight 3PN Einsteinian 
values. The 10PN data are particularly interesting, con- 
sidering that they may be regarded as good approximants 



= 0.1371 



{{E GR )2PN 

LSO 



If the space transformation which relates the Brill-Lindquist so- 
lution with the Misner-Lindquist one through 3PN order l.'j-ll 
is regarded as induced by canonical transformation, the Brill- 
Lindquist and Misner-Lindquist solutions are physically identical 
at this level. 

This determination is not uniq ue f or approximate expressions. 
The method presented in paper |42| is also valuable. It will not 
be employed here since we are interested in the comparison with 
other works. 



(2PN) 



□ 



Skeleton model A: corot. 3PN 
Kidder et al. □: irrot. Einstein 

test particle v — 0: Wilson 
EOB + Pade X : Cook 

Grandclement et al A: Baumgarte 



0.18 

„3/2 



FIG. 1: Position of the last stable circular orbit in the fre- 
quency/energy plane for equal- mass binaries in different mod- 
els. The sequence of black circles corresponds to different 
post-Newtonian truncations of the skeleton energy (from the 
2PN to the 10PN order when moving to the left). The values 
of the 3PN Einsteinian energy for non-spinning black holes are 
computed in Ref. |3U] based on the works of Damour, Jara- 
nowski, Schafer HMlHEl and Blanchet, Faye miMlIJ 
with a; s tatic = 0. The modifications coming from the spins in 
the corotational case are given in Ref. |3Q|. The plus cross 
refers to the point obtained within the hybrid method of Kid- 
der, Will, and Wiseman [4j|, in which terms are added to the 
2PN equations of motion in order to get the exact test-particle 
limit. The data corresponding to the white circle have been 
calculated from a sequence of quasi-equilibrium configurations 
for corotating Misner black holes by Grandclement and col- 
laborators 00] using a Wilson-Mathews-type approximation 
(conformal thin-sandwich approach H). The time cross indi- 
cates the values computed by Cook |25j| for two Misner black 
holes assuming a Bowen- York-like extrinsic curvature (con- 
formal transverse-traceless approach 4j). The white triangle 
shows the position of the last stable circular orbit determined 
by Baumgarte [26| with the help of the puncture method, 
which means in particular that the black holes are of Brill- 
Lindquist type, and that the extrinsic curvature is the one of 
Bowen- York (conformal transverse-traceless approach). The 
star refers to the effective-one-body approach with Pade ap- 
proximant [2^1 . Finally, we have indicated the frequency and 
the energy in the test-body limit {y = 0). 



to the full skeleton model. Remarkably, their position in 
the (ir 3 / 2 , E) plane relative to the 3PN data is nearly the 
same as that of the point referring to numerical simula- 
tions relative to the point referring to the 3PN data in 
general relativity. 

In Fig. ^ the relative location of the last stable orbit 
indicated by the white circle obtained from numerical 
computations with respect to the 3PN result for coro- 
tating bodies (black triangle) is somehow plausible if we 
believe that it is essentially due to the Wilson-Mathews 
truncation. It is indeed comparable to the relative loca- 
tion of irrotational 2PN Einsteinian and Wilson-Mathews 
data. We also want to point out that a realistic estima- 
tion of the last stable orbit does not necessarily follow 
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TABLE I: The last stable circular orbit parameters in the skeleton model truncated at various PN orders for equal-mass 
binaries. 





1PN 


2PN 


3PN 


4PN 


5PN 


6PN 


7PN 


8PN 


9PN 


10PN 


E 


0.5224 
-0.1622 


0.1077 
-0.0689 


0.0725 
-0.0566 


0.0624 
-0.0529 


0.0583 
-0.0514 


0.0563 
-0.0508 


0.0554 
-0.0505 


0.0548 
-0.0503 


0.0546 
-0.0503 


0.0544 
-0.0502 



from the standard 3PN calculation. It may well be close 
to the star-like point derived from the effective one-body 
model with Pade approximant. There is currently no 
definite argument to decide what is the most likely state- 
ment. As a result, we regard the position of these points 
as defining the present uncertainty in the location of the 
last stable orbit derived analytically. 
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APPENDIX A: POLYNOMIAL EQUATIONS FOR 
THE COEFFICIENTS * A IN d = 3 DIMENSIONS 

Equation l|5.2|) for A — 2, after multiplication by 
can be written in the form ^2^1 = A (7) + ^i^i(4)> 



where P 1{7) = %\ + Gp lt Vu/(2r 12 c 3 ) and P 1(4) = 
[Grrii/ (2ri2C 2 )] 2 [^f + p\j '(m 2 c 2 )] are polynomials in *i 
with non-zero constant terms, of order 7 and 4, respec- 
tively. A similar relation holds for A = 1. We now sub- 
tract P2 (7) from the left and right-hand side, and square 
the new identity. We obtain: 

Next, we multiply the equality IjAljl by ( V I / 1 ) 14 = 
*? 8 , make the substitution * 2 ^^(Pim + 

^ffP^l), and expand the terms (P 1{7) + ^fP^)™, 

1 12 

n 6 N* , in powers of P 1 with the help of the 
binomial formula. If n is even, the contribution 

hA n ) — V^™/2 (n\p2k . T/ 8(n/2-fe) p n/2-k .,, f n \ _ 
M l = <Lfc=0 bfcPl(7) W l M(4) ' Wltil \2k) ~ 

n\/[(2k)\(n — 2k)l], is of polynomial kind, whereas that 

, v n/2-l / „ \ p 2fe+lvr.8(n/2-fc-l)+4 p n/2-fc-i+l/2 
01 <Lfc=0 \2k+l) r {7) W l *1(4) na& 

(n) 

the form of a polynomial N± multiplied by the al- 
gebraic function Pyuy M± and have degrees 
max <fe<„/2[7 x 2k + 12 x (n/2 — k)] = 7n and 
max </ c <„/ 2 _i[7x (2fc+l) + 12x (n/2-fc-l)+4] = 7n-3, 
respectively. The same decomposition holds for odd n 
with different expressions for M\ and N\. At last, we 

gather all polynomial quantities into the first member of 

1/2 

the equality while the non polynomial contributions P 1 
are recast to the second member. The final equation is 



7 3 2 

(Ml n2 A 4> \ p2fe ^8(7-/5) p 7-fe Gp 2 iV 2i , nv& 49 \^ f 7 ^ p2fc+1^8(3-fc) p3 -fe , f Gp 2 jV2i \ ^9 

k— k— 



f Gm 2 *>r 2 
V 2ri 2 c 2 



E 

Lfc=0 



12 \ p2k , T ,8(6-fe) p6-k , P2* 



2A- 



nZK ,T/°V"— *) JDO-K 1 
M(7)*l ^1(4) + 



1 (4) T o c 2 2^ 



k=t) 



2k 



p2k , T .8(4-fe) p 4-fc 
r l(7) W l P 



= ~P 



1/2 



1(4) 



V 2r 12 c 2 



1(4) 
3 



(Of-, -n 2 Y^I 14 1 p2fc+l i& 4(13-2A i ) p 6-ft Gp 2i V 2 i ( . 49 ^ / 7 \ 2fe 4(7-2fe) p3 

C* 1 " 1 ) Z^fc + iJ^w *i ^w--^-^ 1 " 1 ^ j^UJ 1(7) 1 p i 

2 



-,4(7-2fc) p3-k 
L(4) 



E 

. fc=0 



12 



2fe+l, T , 4 ( 11 - 2fe )p5-fc _,_ gjjj 8 



1(4) 



m?,c 2 



E 

fc=0 



2fc + l/ jPl ( 7 ) 



2fe+l v j / 4(7-2fe) p 3-fe 



P, 



1(4) 



(A2) 
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We see that the left-hand side of Eq. (| A2|> is made of a sum of polynomials of degrees 100, 99, 98, 98, 98, the first 

1/2 

of them involving the only non-zero constant term. The right-hand side is the product of P 1 ^ and polynomials of 
degrees 97, 96, 95, 95 respectively that cancel when =0. As a consequence, is the root of a polynomial of 
degree 200 with a non-zero constant term. 



APPENDIX B: PROOF OF EQUALITY OF ADM 
AND KOMAR MASSES IN THE CASE OF TWO 
BODIES ALONG CIRCULAR ORBITS 

In this appendix, we restrict ourselves to a binary sys- 
tem. The black holes are thus labeled by A, B, . . . £ 
{1,2}. Our starting point is the system of two sets of 
equations satisfied by the monopoles aA and Pa enter- 
ing the conformal factor \1/ and the auxiliary function 
X = N^: 

PAiVAi 



aA 



3d-2 



Pa 



/ 3d- 2 PAi V Ai 
Pa=XA \ d-2 c * A 



(Bla) 



daA(^A,r 12 J) 
dr 12 



d-2( * A 
a A - Pa 

2ri2 V XA 



(B5b) 



Note that, writing down Eqs. (|B5b|) . we have already 
employed the relations jB5a}. 

The ADM Hamiltonian, restricted to circular orbits, 
reads 



H(ri2,j) = c 2 ^2 a A(ri2,j)- 



(B6) 



A=l 



The circularity condition implies dH(ri2,j)/dr±2 = 0, 
which is equivalent to 



E 

A=l 



da A(r 12, j) 
dr 12 



(B7) 



After substituting l|B6|) into ilB7p we obtain the equality 
(valid along circular orbits): 



Pi 
m A c 2 



E( da A 



A=l 



\8^a dri2 



OT12J 



[1+ d-2m A c** A 



with (here B ^ A) 



*4 = 1 + A t 



XA 



I -A, 



r d-2 
1 12 

Pb 

r d-2 
'12 



(Bib) 



(B2a) 



(B2b) 



It is straightforward to show that for a binary system 
described by the functions W and % of the purely monopo- 
lar form given by Eqs. (|3.14|l and l|4.8(l . the ADM mass 
coincides with its Komar mass if and only if the quanti- 
ties aA and Pa satisfy the simple relation: 



2 

E 

A=l 



(a A -Pa)=0. 



(B3) 



We shall assume that the relative orbits are circular. 
After introducing the reduced relative momentum pi and 
eliminating its square through the relation p 2 — G 2 {mi + 
WI2) 2 "j 2 1 V 2 2 , all quantities become functions of ry2 and j 
only. By virtue of Eqs. (|Bla(l . the effective mass a a can 
then be treated as a function of the three variables *&a, 
r i2 , and j: 

aA^aA{^A(ri2,j),ri2,j). (B4) 
With the help of Eqs. (|B1I) . it is not difficult to check that 
the following relations, crucial for our proof, are fulfilled 
by a A : 

da A (^A,ri2,j) _ Pa 

d^A Xa ' 



(B5a) 



According to Eq. I|B5J) . it can be rewritten as 



A P A d* A d-2 ( * A 

> — « — + "7, PA ola 

^ \XA or 12 2r 12 \ XA , 



= . (B8) 



On the other hand, the equations (|B2a|) in the circular 
case take the form 



*^(ri2,j) = l + A d 



aB(^B(ri2,j),ri2,j) 

„d-2 



B^A. 



By differentiating both sides with respect to ri2 and re- 
sorting to relations (|B5|) . we obtain the system of two 
equations (where B ^ A) 



d^A 

dr 12 



d-2 / * B 
d _i Ad I a B + Pb 



2r 



12 



XB 



A d p B d^ B 
r^ 2 XB dr 12 ' 



which we solve with respect to the derivatives a/ dr\2- 
After eliminating ty A and xa by means of relations (|B2|I , 
the solution reduces to 

^ = -^A d (a B +p B ), B + A. (B9) 
ori2 2rf 2 



To complete the proof of Eq. I|B3(I . it is enough to 
substitute Eqs. (|B9(I into the left-hand side of equality 
and make use of Eqs. (|B2(1 . 
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APPENDIX C: ENERGY ALONG CIRULAR ORBITS AS A FUNCTION OF THE PARAMETER x UP 

TO 10PN ORDER 

The coefficients ei (for i = 5, . . . , 11) entering formula (|6.1()(l for energy E along cirular orbits as a function of the 
parameter x read 



3969 105553 799673 2 502145 , 135247 

e*. = 1 v v H v H 

256 768 2304 3456 62208 

45927 662619 233737 , 1162513 o 56481 4 3777 

e 6 = 1 v v -\ v A H v 

1024 1024 128 1024 1024 1024 

264627 

G7 — 

2048 



16343855 136513223 , 67182599 , 249565613 4 766814851 * 148888223 R 

v v 1 H v A H z/ 4 lA 

6144 18432 27648 31104 497664 13436928 

12196899 993480761 403684801 2 19713422131 3 532659868169 4 1096304004055 5 
32768 + 98304 V 16384 V 884736 V + 3981312 " 11943936 V 

42034682027 6 29197027769 7 
+ 17915904 V + 644972544 V ' 

70366725 2394150375 13003419935 2 9560242585 3 248408863105 4 129128080475 5 
65536 + 65536 ^ 196608 ^ 32768 ^ + 196608 V 98304 ^ 

30481497595 fi 189766295 7 1599565 8 

H ^ b H v 7 H j 

98304 16384 65536 

813439341 100457157571 69694465361 2 44773413415657 3 558727148465215 
262144 + 786432 V 589824 V 21233664 V + 63700992 1 



' "' 262144 + 786432 V 589824 ^ 21233664 V + 63700992 U 

2097054471026245 5 500159286898463 6 858674754130903 7 233795515077109 
191102976 V + 143327232 V + 2579890176 V 5159780352 

177721652120353 9 
417942208512 V ' 

4710988269 228788885795 411952929451 2 28395688669829 3 3163918218071455 4 
611 ~ 524288 + 524288 V + 4718592 V 2359296 V + 63700992 V 

8137959523310875 5 6129421057080619 6 9328338011915891 7 180325112745666065 
127401984 V + 1146617856 V + 322486272 V 30958682112 V 

11728912802084585 9 1098645689138995 10 
278628139008 V + 2507653251072 V ' 
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